// Model 

var 
    c, y, g, k, n, w, R, lam, phi_sh, b_sh, pinf, i, mc, mu_sh 
	lam_f, R_f, c_f, n_f, w_f, y_f, k_f 
	iobs, pinfobs, dy, dc, dy_f, dinve, inve 
	e i_sh bn rir ;

varexo 
    eps_phi, eps_b, eps_t, eps_i, eps_g, eps_mu ;

parameters 
    rho_phi, rho_b, rho_th, rho_g, rho_mu,
    sig_phi, sig_b, sig_t, sig_i, sig_g, sig_mu 
    slop_p psi_i, psi_p, psi_g, psi_y z 
    phi_p, sig, psi, bet, alph, del 
    theta, PISS, MC, phi_k 
    phi, v, b, a, tau, lgv ;

sig   = 2 ;
psi   = 2 ;
bet   = 0.9935^(1/4) ; 
alph  = 1/3 ;
del   = 1-(1-0.08)^(1/4) ; 
theta = 5 ;
PISS  = 1.005 ; 
phi_k = 40 ;

load trend_params

phi = trend_params(1) ;
v   = trend_params(2) ;
b   = trend_params(3) ;
a   = trend_params(4) ;
tau = trend_params(5) ;
lgv = trend_params(6) ;

psi_i = 0.86 ; // JPT (2011)
psi_p = 1.71 ; // JPT (2011)
psi_g = 0.21 ; // JPT (2011)
psi_y = 0.05 ; // JPT (2011)

load est_params

rho_phi = est_params(1) ;
rho_b   = est_params(2) ;
rho_th  = est_params(3) ;
rho_g   = est_params(4) ;
rho_mu  = est_params(5) ;
sig_phi = est_params(6) ;
sig_b   = est_params(7) ;
sig_t   = est_params(8) ;
sig_i   = est_params(9) ;
sig_g   = est_params(10) ;
sig_mu  = est_params(11) ;

z       = 1.0036 ;
slop_p  = 0.01 ;

iSS   = z*PISS*lgv/bet ;
RSS   = z*lgv/bet;
MC    = (theta-1)/theta ;
NoK   = ((RSS-1+del)*(1/MC * (1/alph) * (a*b*z)^(alph-1)))^(1/(1-alph)) ;
WSS   = MC*(1-alph)*a^(1-alph)*(z*b)^(-alph)*(NoK)^(-alph) ;
YoK   = (RSS-1+del)*(1/MC*(1/alph)*(1/z)) ;
CoK   = (1-0.2)*YoK-1+(1-del)/z ;
KSS   = (WSS*(1-tau)*phi/(v*b^psi)*(NoK)^(-psi)*(CoK)^(-sig))^(1/(psi+sig)) ;
NSS   = KSS*NoK ;
CSS   = KSS*CoK ;
YSS   = KSS*YoK ;
lamSS = phi*CSS^(-psi)/lgv ;
GSS   = log(0.2*YSS) ;
inveSS = KSS*(1 - (1-del)/z) ;
NSS   = KSS*NoK ;

phi_p = (theta-1) / slop_p ;

model ;

    // Sticky price economy
    
    (i/(z*PISS*lgv/bet)) = (i(-1)/(z*PISS*lgv/bet))^psi_i * (pinf/(PISS))^(1*(1-psi_i)*psi_p) * (exp(dy)/exp(dy_f))^(1*psi_g) * ((y)/(y_f))^((1-psi_i)*psi_y) * exp(i_sh) ;

	lgv / bet * (1 + phi_k*exp(mu_sh)*((k/k(-1)*z)*exp(mu_sh)-1))*(phi_k/2*(z^2*(k(+1)^2 / k^2)*exp(mu_sh)^2 - 1)+R(+1))^(-1) = (1/z) * lam(+1) / lam ; 
    exp(phi_sh) * phi * c^(-sig) = lam ;
    exp(phi_sh) * v * (b * n)^psi = w * lam * (1-tau) ;
    y = exp(b_sh)^(1-alph)*(a)^(1-alph) * (k(-1)/z)^alph * (b*n)^(1-alph) ;
    y = c + k + exp(g) - (1-del)*k(-1)/z + (phi_p/2) * y * (pinf/pinf(-1) - 1)^2 + (phi_k/2) * (1/z) * k(-1) * (k / k(-1)*z - 1)^2 ;
    w = mc * (1-alph) * y/(b*n) ;
    R = mc * alph * y/k(-1) * z + 1 - del ;

    lgv * pinf(+1) * lam = lam(+1) * i * bet * (1/z) ;
    log(pinf/PISS) = bet / lgv * log(pinf(+1)/PISS) + slop_p * (log(mc/MC)) - e ;

    phi_sh = rho_phi * phi_sh(-1) + (1/100)*sig_phi*eps_phi ;
    b_sh   = rho_b * b_sh(-1) + (1/100)*sig_b*eps_b ;
    e      = rho_th * e(-1) + (1/100)*sig_t*eps_t ;
    i_sh   = (1/100)*sig_i*eps_i ;
    log(exp(g)/y) = (1-rho_g)*(log(0.2)) + rho_g * log(exp(g(-1))/y(-1)) + (1/100)*sig_g*eps_g ;
    mu_sh  = rho_mu * mu_sh(-1) + (1/100)*sig_mu*eps_mu ;

    // Flexible price economy

	lgv / bet * (1 + exp(mu_sh)*phi_k*(z*(k_f/k_f(-1))*exp(mu_sh)-1))*(1*phi_k/2*(z^2*(k_f(+1)^2 / k_f^2)*exp(mu_sh)^2 - 1)+R_f(+1))^(-1) = (1/z) * lam_f(+1) / lam_f ; 
    exp(phi_sh) * phi * c_f^(-sig) = lam_f ;
    exp(phi_sh) * v * (b * n_f)^psi = w_f * lam_f * (1-tau) ;
    y_f = exp(b_sh)^(1-alph) * (a)^(1-alph) * (k_f(-1)/z)^alph * (b * n_f)^(1-alph) ;
    y_f = c_f + k_f + exp(g) - (1-del)*k_f(-1)/z + 1*(phi_k/2) * k_f(-1)/z * (z*k_f / k_f(-1) - 1)^2 ;
    w_f = (theta-1)/theta * (1-alph) * y_f/(b*n_f) ;
    R_f = (theta-1)/theta * alph * y_f/k_f(-1) * z + 1 - del ;

    // Observation equations

    dy      = 1*log(z*y/y(-1)) ;
    dc      = 1*log(z*c/c(-1)) ;
    dinve   = 1*log(z*inve/inve(-1)) ;
    iobs    = 1*log(i/1) ;
    pinfobs = 1*log(pinf) ;

    // Additional equations
    inve = k - (1-del)*k(-1)/z + (phi_k/2) * (1/z) * k(-1) * (k / k(-1) - 1)^2 ;
    dy_f = log(z*y_f/y_f(-1)) ;
    bn   = (b*n) ;
    rir  = log(i)-log(pinf(+1)) ;

end ;

initval ;
    c       = CSS ;
    y       = YSS ;
    g       = GSS ;
    k       = KSS ;
    n       = NSS ;
    w       = WSS ;
    R       = RSS ;
    lam     = lamSS ;
    phi_sh  = 0 ;
    b_sh    = 0 ;
    pinf    = PISS ;
    i 	    = iSS ;
    mc 	    = MC ;
    mu_sh   = 0 ;
    lam_f   = lamSS ;
    R_f     = RSS ;
    c_f     = CSS ;
    n_f		= NSS ;
    w_f		= WSS ;
    y_f		= YSS ;
    k_f 	= KSS ;
    iobs	= log(iSS) ;
    pinfobs = log(PISS) ;
    dy		= log(z) ;
    dc		= log(z) ;
    dy_f	= log(z) ;
    dinve	= log(z) ;
    inve 	= inveSS ;
    e  		= 0 ;
    i_sh  	= 0 ;
    bn 		= b*NSS ;
    rir     = log(iSS)-log(PISS) ;
end ;

steady(maxit=30000,solve_algo=4) ;

shocks ;
    var eps_phi = 1^2 ;
    var eps_b   = 1^2 ;
    var eps_t   = 1^2 ;
    var eps_i   = 1^2 ;
    var eps_g   = 1^2 ;
    var eps_mu  = 1^2 ;
end ;

stoch_simul(order=1, periods=0, nograph, 
    conditional_variance_decomposition=[2, 4]) ;
